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Among the many interesting features displayed by graphene, one of the most attractive is the 
simplicity with which its electronic structure can be described. The study of its physical properties 
is significantly simplified by the linear dispersion relation of electrons in a narrow range around 
the Fermi level. Unfortunately, the mathematical simplicity of graphene electrons is only limited 
to this narrow energy region and is not very practical when dealing with problems that involve 
energies outside the linear dispersion part of the spectrum. In this communication we remedy this 
limitation by deriving a set of closed-form analytical expressions for the real-space single-electron 
Green function of graphene which is valid across a large fraction of the energy spectrum. By 
extending to a wider energy range the simplicity with which graphene electrons are described, it is 
now possible to derive more mathematically transparent and insightful expressions for a number of 
physical properties that involve higher energy scales. The power of this new formalism is illustrated 
in the case of the magnetic (RKKY) interaction in graphene. 

PACS numbers: 



I. INTRODUCTION 



Graphene-related materials have been in the scientific 
limelight for the past few years due to the numerous ap- 
plications envisaged for them[l|. Besides the huge poten- 
tial for applicability, one key feature that makes graphene 
particularly popular is the simplicity with which many of 
its physical properties can be described, primarily due to 
the simple dispersion relation for its electrons. The lin- 
earity of this dispersion relation around the Fermi level 
enables the description of graphene electrons in terms of 
massless Dirac fermions 2]. This introduces a great level 
of mathematical transparency in the portrayal of their 
properties, which nevertheless is limited only to a narrow 
range of energies around the Fermi level. Energy values 
outside this range are often needed, for example when 
gated graphene systems are considered [l[ or when calcu- 
lation of a relevant physical quantity requires an integral 
over energy, but lack the mathematical transparency of 
those within the linear dispersion regime. 

In this communication we show how this limitation 
can be circumvented by deriving a fully analytical closed- 
form expression for the single-electron Green function of 
graphene in real space, a quantity that is instrumental 
in describing the behaviour of graphene electrons. Be- 
cause Green functions are used in the study of several 
physical properties, improvements in their mathematical 
description will enable far more transparent and insight- 
ful expressions for the corresponding physical quantities. 
This is particularly true for distant-dependent interac- 
tions across a graphene sheet, for example the effect of 
an impurity on the physical properties of the system at 
a certain distance away from where it is introduced. An- 
other example is the interaction of two impurities embed- 



ded into the sheet. In both cases, for distances of more 
than a few lattice spacings, the interactions involved are 
mediated by the conduction electrons of the graphene 
host. The Green functions calculated in this paper de- 
scribe the equilibrium properties of these electrons at low 
temperatures, allowing us to investigate the underlying 
interactions within a mathematically transparent frame- 
work. Such a methodology allows for the prediction of 
certain features of the interaction without recourse to 
numerical calculations. 

The remainder of the paper is organised as follows. 
The general method for calculating the Green function 
required is introduced in Section |TTJ The important di- 
rections of the graphene geometry, namely the armchair 
and zigzag directions are illustrated and an explicit cal- 
culation for the real-space off-diagonal Green function el- 
ement in each of these directions is performed in sections 
III Al and III Bl respectively. The accuracy of our approach 
is demonstrated by comparison with a fully numerical 
calculation. An extension of the method for arbitrary 
directions and inter-sublattice cases is discussed in sec- 
tion III CI The potential applications of our approach are 
discussed in Section lllll before an explicit illustration is 
given for the case of the magnetic interaction in graphene 
in section IIII Al Here a fully analytical method is used 
to derive the principal features of the interaction within 
the RKKY approximation. 



II. METHOD AND CALCULATIONS 

The general formula for the single-electron Green func- 
tion, within the nearest-neighbour tight-binding frame- 
work, in its eigenstate basis is 
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FIG. 1: Schematic of the graphene lattice with the zigzag 
(armchair) direction denoted by the blue (red) arrow marked 
'Z' ('A'). The filled and hollow circles represent carbon atoms 
on each of the two sublattices which compose the graphene 
lattice. 



where E is the energy, \k, ±) is the eigenvector labeled by 
the wave vector k and £ ± is the corresponding eigenvalue 
defined as 



£±(k) = ±iWl + 4cos( 



^^ ya )cos(k x ^) +4cos 2 (fc 3: ^) . 

(2) 

The quantities a and t correspond to the lattice pa- 
rameter of graphene and its nearest-neighbour electronic 
hopping Q, respectively, which are hereafter used as our 
units of distance and energy. Our choice is such that 
the x (y) direction is aligned to the zigzag (armchair) 
geometry of the graphene lattice as is shown schemati- 
cally in Fig. [TJ A simple unitary transformation defines 
the real-space basis | j, C) , where the index j labels the 
two-atom unit cell and the index £ refers to the intra- 
cell atoms corresponding to the two distinct sublattices 
of graphene, represented by filled and hollow circles in 
Figure [T] When projected onto two different states \j, () 
and \f , C) located in real space by the respective vectors 
Rj and Rji , the Green function is written as 
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where the integrals over k x and k y are performed over 
the first Brillouin Zone of graphene. Although we have 
selected two states that belong to the same sublattice 
(C = £'), this constraint can be easily relaxed and gener- 
alized to describe the propagator between sites on dif- 
ferent sublattices. Before proceeding, we will outline 
the basic steps taken to obtain the Green function. We 
tackle the first integral by analytical continuation to the 
complex plane, where it is subsequently solved using the 
residue theorem. 0, [E| The remaining integral can then be 
solved in the case of moderately large separation vectors 
by using the Stationary Phase Approximation (SPA). [6] 



A. Armchair Direction 

We first consider the case of separation vectors Rj—Ry 
along the (armchair) y-direction. By showing how to ob- 
tain the Green function for this particular case we hope to 
illustrate the general method for calculating Green func- 
tions for any direction. Note that the position vectors 
appear only as a difference and can be further simplified 
by defining A = \Rj — Rj>\- In this case, the integral is 
performed over the Brillouin Zone shown in Fig. [5] and 
the first integral, over k y , can be evaluated by extending 
k y to the realm of complex numbers and changing the in- 
tegration contour from a straight line on the real axis to 
the boundaries of a semi-infinite rectangle in the upper 
half of the complex plane, with its base lying on the real 
axis between — ^ and . Because the integrand van- 

av3 av3 

ishes in the limit Im[fcj,] — ¥ oo and because the parts of 
the contour that are parallel to the imaginary axis cancel 
each other out, the fc^-integral can be evaluated by sim- 
ply identifying the poles of the integrand lying inside the 
integration contour and finding their respective residues 
3, that is, 
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Note that the scalar product (j, (\G{E)\j' ,('} is now more 
concisely expressed as Qa (E) and that 
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defines the wave vector q that comes out of the first in- 
tegral. Although Eq.© provides a unique definition for 

cos ( ^ ), it does not specify the sign of q uniquely. Its 
sign is selected by imposing that q must necessarily lie 
within the integration contour of the fc y -integral. 

The fc^-dependence contained in the wave vector q, as 
seen in Eq.©, means that the integrand in Eq. (|4|) is 
an oscillatory function of k x that oscillates very rapidly 
for large values of the separation A. In this case, the 
only non-vanishing contribution to the Green function 
comes from regions for which the phase of the exponential 
function is stationary. To locate these stationary points 
we must impose that dq/dk x = 0, which leads to solutions 
of the form 
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Note that due to a topological change in the constant 
energy surfaces of the function £+{k) at E = ±t we sep- 
arate the energy band into two separate regions, namely, 
the inner region defined by \E\ < \t\ and the outer region 
defined by \E\ > \t\. Both stationary solutions written 
above are valid throughout the entire energy spectrum. 
However, outside those specified regions, albeit solutions, 
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FIG. 2: Constant energy plots of the function £+(fc) in recip- 
rocal space for a few different energies. Horizontal and vertical 
axes are rescaled as k x a/2 and k y ay/3/2, respectively, so that 
they are plotted as dimensionless quantities. The rectangular 
shaded area delimits the first Brillouin zone over which the 
integrals for the armchair direction case are taken. Constant 
energy plots for two specific energies, E — 0.7\t\ (solid) and 
E — l.8\t\ (dashed), are drawn with thicker lines. The cor- 
responding stationary wave vectors q for these energies are 
highlighted with arrows. Note that the sign of q, shown as 
positive here for simplicity, must be chosen according to the 
conditions outlined in the text. 



they give rise to complex values for the wave vectors 
q. With complex wave vectors, the integrand of Eq.Q 
tends to vanish for any sizable separation A, meaning 
that the stationary values outside the ranges listed in 
Eq. © should have very little influence in the overall re- 
sults for the Green functions. Q The same can be under- 
stood from purely geometrical arguments applied to the 
constant energy surfaces of £+{k) in reciprocal space, de- 
picted in Fig. [3J When searching for stationary solutions 
for q, which in this case lie parallel to the y-axis, the 
two solutions resulting from Eq.© are the only possible 
(real) ones within the rectangular Brillouin zone of the 
hexagonal lattice. The tilde symbol (~) will hereafter 
be used to refer to the values of k x and q satisfying the 
stationary condition. Therefore, the wave vector q when 
expanded in a Taylor series around the stationary value 
k x has no linear component and, up to second order, is 
approximated by 
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Note that the sign of Cx must be chosen to ensure that 
q lies within the k y integration contour as before. The 
sign of C2 is determined by its correspondence to the 
curvature of q at k x . If we now insert Eq.© into Eq.Q 
and make use of the fact that k x and q will not vary very 
much around their respective stationary values k x and 
q, we are left with a much simplified expression for the 
Green function, which now reads 
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The remaining integral is a well known Gaussian integral 
whose solution gives 
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This can be rewritten in a more transparent fashion 
using the definitions in Eqs. (O, © and (|9|) to provide 
a completely analytical expression for the off-diagonal 
Green function matrix element between two graphene 
sites separated by a distance A along the armchair di- 
rection. For positive values of energy (E > 0), we find 



Qa(E) = 



(E 2 + 3t 2 ) Vt 2 - E 2 



iVE( iE +^f^ ) A ' 

E ,E 2 -5t 



y/VE 2 -9t 2 



l ^f(t 2 -E 2 ){E 2 -9t 2 ) ,^ 
it 2 > 



if E < \t\ 
if E > Itl 



(12) 



where A' = -^ A =. We again consider the distinct cases 

av3 

E < \t\ and E > \t\ and note that the only occa- 
sion when both stationary points contribute is at en- 
ergies very near E = ±t [7j. Fig. [3] compares both 



the real and imaginary parts of the Green function for 
the case of A = lOv^a obtained by the analytical ex- 
pression above with those obtained through a numeri- 
cal evaluation of Eq. (J3j . For E < we note that the 
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FIG. 3: Qa(E) as a function of energy (in units of \t\) for 
the case of A = Wy/Sa. Top (Bottom) panel shows the real 
(imaginary) part of the Green function. Lines correspond to 
the results evaluated by Eq. (|ll[) whereas points are the result 
of brute-force numerical calculations of Eq.©. 
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FIG. 4: Fraction T\% of the bandwidth for which the error 
between the analytical and brute force numerical results is 
below 1%, plotted as a function of the separation A (in units 
of a) . The circles (triangles) correspond to separations in the 
armchair (zigzag) direction. 



real (imaginary) part of the Green function is an odd 
(even) function of energy, as can be seen in Fig. [3J and 
make use of the relations Re(Q A (— E)) = —He(Q / \(E)) 
and Im(g A (-E)) = Im(g A (E)). 

At first glance one might think that the replacement 
of the well established linear dispersion approximation 
with another that is valid only for asymptotically large 
values of separation is unlikely to improve the range of 
validity of the overall result. However, as seen in Fig. 
131 there is hardly any noticeable difference between the 
analytical and numerical results across the entire energy 
band. Because the analytical expression relies on the 
stationary phase argument, this remarkable agreement is 
likely to regress as the separation (A) is reduced. To 
test how good an approximation Eq. fTTl) is, in Fig. H]the 
fraction T\% of the bandwidth for which the relative er- 
ror between the numerical and analytical evaluations is 
less than 1% is plotted as a function of the separation. 
The plot with circular points corresponds to the arm- 



chair direction. Even for small separations (A ss 10a), 
the energy range for which the Green functions are very 
accurately described by our analytical expression exceeds 
90% of the bandwidth. In other words, there is only 
a very narrow energy range in which the disagreement 
exceeds 1%. Most remarkably, as the separation is in- 
creased this small range decreases very rapidly, indicat- 
ing that Eq. (ITT1) is capable of accurately describing the 
Green function across almost the entire energy spectrum 
for separation values larger than a few lattice parame- 
ters. This is a major advantage when compared to the 
narrow fraction of the bandwidth that meets the linear 
dispersion criterion. 



B. Zigzag direction 

We now turn our attention to the case of separation 
vectors Rj — Rj> along the (zigzag) cc-direction. By fol- 
lowing the same approach described for the armchair di- 
rection the relevant Green function can be similarly cal- 
culated. The major difference between the calculations 
is that the ordering of the integrals is swapped. For the 
zigzag direction, we first perform a contour integral over 
k x before making use of the SPA to solve the remaining 
integral over k y . We make a different choice of Brillouin 
Zone, as shown in Fig. [SJ which will simplify the selec- 
tion of stationary points later. By performing the first 
integral similarly to before, we arrive at an expression 
analogous to Eq. for the off-diagonal Green function 
in the zigzag direction 
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where q now represents the poles coming out of the k x 
integral, which are given by 
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It should be noted that in the zigzag case there are two 
contributions to the integral arising from the two possible 
sign choices of the poles in the definition above. The 
correct overall sign for q in each case is selected as before 
by ensuring that q lies within the integration contour. 
The contributions from each pole must then be summed 
to give the final result. As before, we assert that the 
only non- vanishing contributions to the integral in Eq. 
(fl~3| occur when the phase of the exponential term is 
stationary. Imposing dq/dk y — we find the stationary 
solution k y = 0. 

Unlike the stationary points calculated for the arm- 
chair direction, the zigzag direction stationary points are 
independent of energy. This fact is evident when the 
constant energy plots in Fig. [S] are examined from the 
perspective of stationary values of k y . The separation of 
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FIG. 5: The constant energy surfaces of the function £+(k) 
as before with the Brillouin Zone for the zigzag case now 
illustrated by the shaded area. The thick lines once more refer 
to constant energy plots for E = 0.7|i| (solid) and E = 1.8|t| 
(dashed). The corresponding stationary wave vectors q for 
these energies are highlighted with arrows. Note that the 
sign of q, shown as positive (solid) or negative (dashed) only 
for simplicity, must be chosen according to the conditions 
outlined in the text. For clarity the arrows are shifted slightly 
in the vertical direction, but it should be noted that all the 
stationary points occur at k y — 0. 



the energy band into two separate regions is not neces- 
sary in this case as the stationary points for both regions 
occur for the same value of k y . The wavevector q is now 
Taylor expanded as before and we find expressions for C\ 
and Ci analogous to Eqs. © and ^ 
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The superscript sign in the expressions for C\ and Ci 
refer to the choice of sign in Eq. [14] Using these ex- 
pressions the integral once more reduces to a Gaussian 
integral whose solution gives 
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where the sum over a includes the contributions arising 
from the choice of sign for the poles. Eqs. (fl5l - IT7|) can 
be combined as before to provide a single fully analytical 
expression for the off-diagonal Green function matrix ele- 
ment between two graphene sites separated by a distance 
A in the zigzag direction. This is found to be 
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where A' = Here we have a single expression that 
describes the Green function across the entire band. In 
Fig. [5] we compare the expression for the Green function 
calculated here with the result of a fully numerical cal- 
culation, for the case of A = 20a. As with the armchair 
direction, an excellent match is found across the entire 
band. The plot in Fig. 2] (triangular symbols) shows the 
discrepancy between the numerical and analytical results 
as a function of distance. Once more we find that beyond 
a couple of lattice spacings there is only a very narrow 
energy range in which the disagreement exceeds 1%. 



C. Other directions and cases 

Having presented the derivation of the Green function 
for the separation vector along the armchair and zigzag 
directions, it is straightforward to generalize it to other 
cases. For arbitrary directions, although the procedure is 
similar, we shall see that the identification of the poles or 
stationary solutions may result from high order polyno- 
mial equations that are not always analytically solvable. 

In the armchair (zigzag) case, the expression for the 
stationary point k x (k y ) is given by an easily solvable ex- 
pression of the form dq/dk = 0. This expression arises 
from the decision to take the contour integral along the k- 
space direction parallel to the separation vector Rj — Rji , 
which results in a phase term equal to the pole of the con- 
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FIG. 6: Qa(E) as a function of energy (in units of \t\) for 
the case of A = 20a. Top (Bottom) panel shows the real 
(imaginary) part of the Green function. Lines correspond to 
the results evaluated by Eq. (|17p whereas points are the result 
of brute-force numerical calculations of Eq.(|3}- 



tour integral. Since the expressions for the poles in the 
armchair and zigzag directions, Eqs. (J5J) and (|14p respec- 
tively, are easily found from Eq. ©, the calculation of 
all the necessary quantities is quite straightforward. To 
extend this approach to an arbitrary separation vector, 
we must first rewrite Eq. ([2]) in terms of fc-space vectors 
fcy and k± which are parallel and perpendicular respec- 
tively to the required separation vector. Following this, 
we must perform the contour integral over fey to get an 
expression for the Green function analogous to Eq. (fj|. 
However, the expression for the poles of this contour in- 
tegral will depend strongly on the separation vector cho- 
sen and will usually result from a high order polynomial 
equation that may need to be solved numerically. It is 
important to note that this equation will only depend 
on the direction, and not the length, of the separation 
vector, so that once the Green function for a particular 
direction has been constructed it is valid for any required 
distance across the graphene lattice in that direction. 

A similar methodology holds for the case of Green 
functions between sites on the different sublattices of 
graphene. In this case Eq. ([3]) must be altered slightly 
to read 
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where now we have £ 7^ with f(k) 

2cos(-^)e 2 ^3 . The integral can now be split into two 
parts with different phase terms coming from the two 
components of f(k). These can then be solved individ- 
ually using the approach described above to give the re- 
quired Green function. 



III. APPLICATION 

The usefulness of having an analytical expression for 
the real space Green function, valid throughout the entire 
electronic energy band, becomes obvious when physical 
properties of graphene involving energy scales outside the 
linear dispersion region are investigated. This is partic- 
ularly advantageous when such properties carry size or 
position dependence because in this case Eq. (ITll for the 
armchair direction, or Eq.(|17p for the zigzag direction, 
can be more concisely expressed in the form 
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so that the E- and A-dependences are clearly distin- 
guished. Furthermore, even in the case when the func- 
tional form of the coefficient A(E) is not particularly sim- 
ple, the expressions in Eqs. (fTTl) and (TTTf can be used to 
expand the Green function in a polynomial series, which 
is undoubtedly far simpler and more treatable than the 
original expression in Eq. © . We anticipate that the abil- 
ity to clearly isolate the distance dependence in the Green 
function will allow a more transparent treatment of some 
of the more eagerly investigated properties of graphene. 

This approach will be shown more clearly in the next 
section when our formalism is used to investigate the 
magnetic interaction between two magnetic moments in 
graphene. This type of interaction is perfectly suited 
for investigation using our approach since it is mediated 
by the conduction electrons of the graphene host. How- 
ever, there are many other interesting physical properties 
that can be explored. The interaction between process- 
ing magnetic moments is one area of particular interest. 
Within the random phase approximation this dynamic 
coupling can be described by an integral over a complex 
function involving the convolution of Green functions. 
Initial numerical results in carbon nanotubes, which are 
closely related to graphene, suggest that the range of the 
dynamic interaction may be greater than that of the more 
familiar static case |19l - [2l| . We anticipate that an exten- 
sion of the description provided below for the static case 
may be useful in attempting to solve the required integral 
analytically and understand the distance dependence of 
the dynamic coupling. The ability of our approach to 
obtain Green functions over a very large fraction of the 
energy band becomes increasingly important in this case 
due to the convolution of Green functions of different en- 
ergies that appears in the expression. 

Another topic that lends itself to our approach is the 
effect of disorder [22J and in particular, the effect of an 
impurity, on the properties of graphene. Friedel os- 
cillations, occuring in the local density of states as a 
function of distance from an introduced impurity, have 
been studied within the linear dispersion regime using a 
Green function formalism (23|. Although the local den- 
sity of states is associated with the diagonal term of the 
Green function, the distance dependence of the oscilla- 
tions will be determined solely by the off-diagonal term 
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calculated here. Similarly, the signatures of magnetic 
adatoms in graphene when probed by scanning tunnel- 
ing spectroscopy have also been investigated using a the- 
oretical approach [2^. This method again avails of Green 
function methods within the linear dispersion regime. We 
anticipate that our approach may allow for an extension 
of such studies to energies beyond the linear dispersion 
regime. 



A. Application to RKKY interaction 

To demonstrate the power and applicability of our 
new formalism, we turn our attention to the magnetic 
interaction in graphene. This interaction determines 
the relative orientation of magnetic moments embedded 
in graphene and has been the subject of many recent 
papers as an understanding of this interaction is 

a major step in the implementation of graphene devices 
in the field of spintronics. When the linear dispersion 
approximation is used, a cut-off function is required to 
prevent the result diverging due to high energy contri- 
butions. There has been some debate about the effect 
of the cutoff function chosen on the resultant interaction 
calculated [8, 9,T2j. Other approaches to circumvent this 
problem involve numerical calculations fl2l [lij which can 
lack the transparency of an analytical solution. Here we 
shall show that the decay rate and oscillation period of 
the interaction as a function of distance emerge natu- 
rally from a simple calculation using our formalism and 
without resorting to an energy cut-off or a restriction 
to the linear dispersion regime. The exchange energy, 
J, within the RKKY approximation [l^| - fl7l ] is propor- 
tional to the static susceptibility, x, which can in turn 
be written in terms of Green functions, allowing us to 
write Ja(Ef) ~ Im J &Ef{E)Q\(E) for two moments 
occupying like-sites in the graphene lattice separated by 
a distance A, where f(E) is the Fermi function. This 
quantity relates to the energy difference between the fer- 
romagnetic and antiferromagnetic alignment of the mo- 
ments, with its sign giving the energetically favourable 
alignment. Using the expression in Eq. (|20[) we write 



around the Fermi energy gives 
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using B^ (Cj ) to denote the Ith derivative of B (Ci) 
evaluated at the Fermi energy. This can be rewritten as 
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In the low temperature limit, T — > 0, this expression 
simplifies to one of the form 

2id(E F )A 

J A (E F )~Y, B{l) ( E ^ Al+ 2 ■ ( 24 ) 

In this form the oscillation period and decay rate of 
the interaction at different Fermi energies can be eas- 
ily extracted. The decay rate in the asymptotic limit [l8[ 
is determined by the leading term in Eq (1241) . namely 
I = 0, suggesting that, in general, J ~ A~ 2 . However, 
at Ep = 0, it is straightforward to determine from Eq. 
(fT2l for the armchair direction, and from Eq. (TTSl for 
the zigzag direction, that the coefficient B^ vanishes 
and the decay rate is in fact determined by the first sur- 
viving term, 1 = 1, resulting in J ~ A -3 for undoped 
graphene, as reported elsewhere H, [IH, Also, at 
Ep = 0, the oscillation period is perfectly commensurate 
with the graphene lattice spacing and thus oscillations 
are masked. When £j? / 0, the leading term does not 
vanish, and the oscillation period is no longer commen- 
surate with the lattice spacing, leading to the observed 
oscillatory interaction [9j that decays as J ~ A~ 2 . Note 
that these conclusions are reached for values of Ep re- 
gardless of whether they lie within the linear dispersion 
regime. The correct decay rate and oscillatory behaviour 
for the RKKY interaction in graphene have emerged nat- 
urally and in a mathematically transparent fashion from 
our formalism, without resorting to the linear response 
approximation or the need for a cut-off function. As far 
as the authors are aware, this is the first time this has 
been performed within a fully analytical framework. 



IV. CONCLUSIONS 



where B(E) = A 2 {E), (3 = T being the temper- 

ature and Ub the Boltzmann constant. The integral in 
Eq (f2~Tj) can be solved by replacing it with a contour 
integral in the energy upper-half plane. In this case 
the poles are given by the zeroes of the denominator of 
the Fermi function, namely the Matsubara frequencies, 
E p = Ep + i(2p + l)irksT where p is an integer labelling 
the poles. Writing the coefficient B(E) as a Taylor se- 
ries, and the wavevector C\(E) as a first order expansion, 



In summary, we have derived closed-form expressions 
for the single-electron Green function of graphene in real 
space that does not rely on the linearity of its disper- 
sion relation near the Fermi level. The full derivation 
of this quantity for the two principal directions investi- 
gated on the graphene lattice has been presented along 
with a discussion for extending the methodology to other 
cases. The expressions are valid across a very large frac- 
tion of the energy band and yet remain mathematically 
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transparent. The newly- acquired simplicity with which 
we can describe the electronic properties of graphene will 
lead to insightful new ways of studying the physical prop- 
erties of this material at energy scales well beyond the 
linear dispersion regime. The approach described here 
for the magnetic interaction can be modified straightfor- 
wardly to extend the validity of expressions derived in 
topics including, but not limited to, modelling the dy- 
namic magnetic coupling and the effects of adatoms in 
graphene systems. 
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